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ABSTRACT 


The  effect  of  round-off  errors  on  the  estimation  of  position  and 
velocity  in  midcourse  navigation  has^seen  analyzed.  The  analysis!  is 
based  on  Kalman's  approach  to  linear  filtering  and  prediction.  -'"Com¬ 
putation  noise""appears  to  be  an  additional  random  force  in  the  dynamic 
system,  and  may  affect  both  convergence  and  equilibrium  of  the  sequen¬ 
tial  estimation  procedure  significantly. 

The  analysis  hili^cn  applied  to  a  satellite  trajectory  estimation 
system.  Axes  and  area  of  the  error  ellipse  (ellipse  of  concentration) 
have- oe an.  expressed  in  terms  of  word  length,  time  interval  between 
observations,  and  number  of  integration  steps. 
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1.  INTRODUCTION 


Selection  of  the  wordlength  in  digital  computers  used  in  guidance 
and  navigation  systems  is  infrequently  considered  to  be  a  serious  design 
problem. 

The  decision  is  most  often  based  on  whether  1  or  2  instructions 
should  be  stored  in  one  memory  word,  and  how  many  digits  of  the  in¬ 
struction  word  should  be  allotted  to  the  order  code  and  address.  Index 
modification  schemes  for  shortening  the  address  may  have  more  effect 
on  the  wordlength  selection  than  accuracy  requirements. 

Reliability  considerations  increase  the  importance  of  wordlength 
studies  significantly,  especially  if  the  digital  computer  is  to  be  used  on 
board  a  space  vehicle  and  operates  independently  of  ground  communica¬ 
tions.  In  optimizing  the  hardware  for  such  a  system,  it  is  necessary  to 
relate  computer  functions  (i.  e.  accuracy)  to  the  performance  of  the  whole 
system,  and  to  design  the  hardware  so  that  intolerable  degradations  in 
performance  re  suiting  from  failures  are  minimized  or  eliminated.  A  typ¬ 
ical  example  of  this  realistic  system  optimization  philosophy  is  the  "Word 
Split  Technique".  It  is  a  technique  whereby  a  digital  word  is  automati¬ 
cally  spiit  after  a  failure  and  the  failing  part  is  excluded  from  the  opera- 
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tion  ’  .  Extensive  studies  at  Hughes  Aircraft  Company  ’  have  shown 

that  application  of  a  "failure  tolerance  concept"  instead  of  a  pure  redun¬ 
dancy  concept  may  very  likely  lead  to  significant  hardware  savings. 
However,  application  of  the  failure  tolerance  concept  to  real  time  data 
requires  a  thorough  analysis  of  the  performance  degradation  caused  by 
dropping  the  least  significant  word  half. 

The  objective  of  this  paper  is  to  analyze  the  effect  of  computer 
accuracy  limitations  on  system  performance  during  midcourse  naviga¬ 
tion.  The  purpose  is  two  fold: 

a.  To  supply  the  tools  for  determining  the  error  contributions 
of  a  (given)  digital  computer  in  midcourse  navigation. 

b.  To  establish  whether  degraded  operational  modes  in  extremely 
long  range  space  missions  may  be  tolerated  , 
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Several  papers  have  been  published  recently  which  describe  the 
application  of  concepts  from  statistical  filter  theory  to  inflight  deter¬ 
mination  of  position  and  velocity  of  a  (manned)  space  vehicle  for  the 
purpose  of  midcourse  guidance.  The  spaceborne  digital  computer  im¬ 
plements  a  dynamic  time  varying  filter  which  weights  the  incoming 
observations  in  an  "optimal"  sense  and  produces  an  up  to  date  optimal 
estimate  of  any  desired  set  of  state  variables. 

The  basic  theory  quite  frequently  applied,  is  described  in  refer¬ 
ences  (1)  and  (2)  where  it  is  shown  that  every  solution  of  the  variance 
equation  converges  to  an  equilibrium  point,  and  that  the  equilibrium 
exists  if  certain  conditions  are  satisfied.  The  variance  equation  is 
described  as  a  stable  computational  method  and  is  expected  to  be  in¬ 
sensitive  to  round-off  errors.  Application  of  the  theorems  and  hypoth¬ 
esis  in  (1)  and  (2)  would  always  lead  to  the  conclusion  that  computation 
errors  are  negligible  compared  to  instrument  anomalies.  It  seems 
questionable,  however,  to  draw  any  conclusions  about  the  propagation 
of  round-off  errors  from  a  theory  which  does  not  permit  such  errors  in 
the  basic  dynamic  model. 

In  the  following,  we  introduce  "computation  noise"  in  addition  to 
the  random  forces  and  measurement  uncertainties,  and  study  its  effect 
on  the  convergence  and  equilibrium  of  the  estimation  procedure. 


2.  SYSTEM  MODEL 


The  fundamental  relations  for  finding  the  best  estimate  of  the 
message  process  in  the  linear  dynamical  system  of  the  form 


CD  x(t+i) 

(2)  Z(t) 
are  given  by 

(3)  &(t+l/t) 

(4)  $::(t+l.t) 

(5)  A*(t) 

(6)  P(t+l,t) 


where 

®(t+l,t) 

x(t) 

Z(t) 

M(t) 

r  (t+i,t) 

v(t) 
w(t) 
x(t/ 1-  1) 


=  $(t+i,t)  x(t)  +r(t+i,t)  w(t) 

=  M(t)  x(t)  +  v(t) 


=  $"'(t+l,t)  x(t/t-l)  +  A*  (t)  Z(t) 

=  $>  (t+  1 ,  t)  -  A‘  (t)  M(t) 

=  ®(t+l,t)  P(t/t-l)  MT(t)[M(t)  P(t/t-l)MT(t)  +  R(t)] 
=  ®(t+l,t)  P(t/t-l)  -  [  P(t/t-l)  MT(t)  ] 

[  M(t)  P(t/t-l)  M(t)  +  R(t)  J  ’ 1  |M(t)  P(t/t-l)  ] 
$T(t+i,t)  +  r(t+i,t)  Q(t)  rT(t+i,t) 


- 1 


=  state  transition  matrix 

=  state  vector 

=  vector  of  measurements 

=  transformation  relating  the  observables  to 
the  state  vector 

=  transition  matrix  for  the  random  force  w(t) 

=  measurement  uncertainty 
=  random  force 

=  estimated  state  vector  based  on  past  observation 
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A  (t) 

P(t+l.t) 

R(t) 

Q  (t) 

It  is  assumed  that 


=  optimum  filter 

=  covariance  matrix  of  the  estimation  error 

=  covariance  matrix  of  the  measurement  uncer¬ 
tainties  v(t) 

=  covariance  matrix  of  the  random  force  w(t) 


(7) 

E[w(t)  ] 

-  E  [v(t)  ]  =  0 

(8) 

E  |w(t)  \vT(t)  ] 

=  Q(t) 

(9) 

E  [v(t)  vT(t)  ] 

=  R(t) 

For  the  derivation  of  the  formulas  see  reference  (1)  and  (2). 
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3.  DEFINITION  OF  COMPUTATION  NOISE 


No  digital  computing  procedure  or  device  can  perform  the 
operations  which  are  its  "elementary"  operations  rigorously  and  fault¬ 
lessly  because  of  the  finite  length  of  a  "digital  word".  Each  time  an 
elementary  operation  is  performed,  a  perturbation  is  introduced  which 
will  cause  a  parameter  to  deviate  from  its  ideal  value.  The  non  ideal 
elementary  operations  will  be  called  "pseudo  operations".  They  form 
a  constantly  renewed  source  of  contamination,  and  their  influence  in¬ 
creases  with  the  number  of  elementary  operations  that  have  to  be  per¬ 
formed.  They  are  therefore  especially  important  in  long  computations, 
involving  many  such  operations.  Long  computations  will  undoubtedly 
be  normal  for  midcourse  estimation  and  decision  making  procedures. 
The  decisive  factor  that  controls  their  effect  is  some  kind  of  stability 
phenomenon.  And  it  is  the  stability  of  the  approximant  procedure  and 
not  of  the  strict  procedure  that  matters.  Estimation  procedures  which 
are  theoretically  identical  (asymptotically) ,  but  differ  in  the  number  and 
sequence  of  pseudo  operations  will  be  characterized  by  a  different  speed 
of  convergence  to  a  different  equilibrium  point. 

Notation  for  pseudo  operations: 


A 

A  +  Tj  where 

hi 

<  2  (K  =  wordlength) 

A©B  - 

A  +  B~+  q 

a 

hu 

|*2-<K+1> 

A.  x  T5  = 

A  .  5  +  T| 

m 

bm| 

j*2-<K+1> 

AfB  = 

A/B  +  qd 

bd ! 

i2-(K+!) 

The  q  's  are  random  variables  and  assumed  to  be  uniformly  dis¬ 
tributed  between  +  2  .  In  other  words,  all  bit  configurations  are 

assumed  to  be  equally  likely  for  -j_2.  k,  where  j  designates  the  truncated 
digits. 
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It  is  remarked  that  scaling  is  also  a  pseudo  operation.  It  is 
identical  to  a  multiplication  or  division,  but  is  carried  out  with  a 
multiplier  or  divisor  outside  of  the  range  of  digital  numbers  in  a  fixed 
point  organization. 

The  pseudo  operations  considered  so  far  are  the  "noisy  operations" 
in  the  failure  free  mode.  These  can  be  characterized  by  a  uniform 
error  distribution  which  converges  rapidly  to  the  normal  distribution 
as  the  number  of  pseudo  operations  performed  in  sequence  increases. 

But  if  permanent  hardware  failures  are  "allowed  to  occur"  in  long 
range  space  missions,  the  computer  may  very  well  degrade  its  per- 
formance^'^’^1^.  The  statistics  of  pseudo  operations  in  degraded 
modes  (time  of  occurrence  is  a  random  variable)  is  not  analyzed  here. 
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4.  EFFECT  OF  COMPUTATION  NOISE  ON  THE  ESTIMATION  ERROR 


The  estimate  of  the  state  vector  is  computed  from 

x(t+l/t)  =  ®D(t+l,t)  X  ^t/t-l)©^  (t)xZ(t) 

where  the  index  D  indicates  that  the  numerical  values  of  the  matrices 
#  .  * 


and  A 

are  erroneous: 

(10) 

>i< 

* 

5 

$ 

+  6  $ 

(11) 

+  6  A 

t> 

a 

n 

A 

The  estimation  error  in  step  (t+l,t)  caused  only  by  computation  noise 
is  therefore  in  a  first  approximation: 

(12)  6x(t+ 1 ,  t)  -  [$xx-$.  xj  +  |  5*  x  -  A*  .  Z  ]  - 

-j  x¥x  x  -A  ,  M  .  x  |  + 

+  fox  +  6A*Z  -  toMx 

The  total  estimation  error  in  step  (t+l,t)  is: 

(13)  xD(t+l/t)  =  x(t+  1/t)  +  6x(t+  1 ,  t) 

For  simplification  it  is  assumed  in  the  following  that  the  errors  intro¬ 
duced  by  pseudo  operations  in  calculating  x(t+l,t)  are  negligible  against 
those  introduced  in  computing  $  and  A  .  Thus, 

(14)  6x(t+ 1 ,  t)  =  64>x  +  6A*Z  -  6A*Mx 


Some  properties  of  6$  and  £A  *  which  follow  from  the  algebra  of 
pseudo-operations  are  listed  below: 


(15)  e[6»] 

|&<I>  x  | 
[dA*x 
[^xK] 
[6A:;:jxK 


E 

E 

E 

E 


=  [  6  A*  ]  =  0 

=  E  [  6$  ]  E  [  x  ]  =  0 


=  0 


=  E  [6$^  ]  e[xK  ] 

=  E  [&AVj  ]  E  xK  ] 

*1 


E  6A 


E  [s$6A*  ]  j-  E  [&$>  ] 

E  |  E  |  6$  J  E  |  J  ,  similar  for  6A' 


Introducing  (2)  and  (13)  into  (12)  we  obtain 


where 


and 


6x(t  + 1 ,  t)  =  6*ft(t/t-l)  +  6A  Mx  (t/t-1) 


x(t/t-l)  =  x(t)  +  x(t/t-  1)  +  6x(t,  t-  1) 


xD(t/t-  1) 


-  x(t/t-l)  +  6x(t,t-l) 


or 


^(t/t-1) 


5A*MxD(t-l/t-2)  +  x(t/t-l)  +  63>x(t-  l/t-2) 
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Repeated  application  of  the  recursion  formulas  above  leads  to: 

(16)  6x(t+ 1 ,  t)  =  2  (6$  +  6A*M)1-1.  6$  x(t-i)  + 

i 

+  £  (6$+  6A*M)X  x(t-i/t-i-l)  + 
i 

+  (6$  +  fiA!,<M)n  6x(t-n) 

The  effect  of  errors  made  in  previous  computation  (observation) 
intervals  is  gradually  forgotten  as  t  — >  oo;  but  the  speed  of  convergence 
depends  on  the  behavior  of  the  state  vector  (first  term  in  16)  and  estima¬ 
tion  error  (second  term  in  16)  as  a  function  of  time.  It  depends  also, 
of  course,  on  the  magnitude  of  6$  and  6A  M.  Assuming  both  to  be  mod¬ 
erately  small,  the  last  term  in  (16)  can  be  neglected. 

Because  of  the  dependency  of  6x  on  the  state  vector,  convergence 
and  equilibrium  are  not  guaranteed  anymore  by  the  conditions  given  in 
[l]and[2]  as  it  will  be  shown  in  a  simple  example  in  section  5. 

It  may  be  sufficient  to  consider  only  first  order  terms: 

(17)  6x(t+l,t)  £  6*x(t)  +  (6S+  6A*M)  x(t/t-l) 
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5.  COVARIANCE  MATRIX  OF  ESTIMATION  ERROR  (VARIANCE  EQUATION 


Applying  the  expectation  operator  in  a  straight  forward  manner 
according  to 

P(t+1)  =  E  [  xD(t+l/t)  x£(t+l/t)  ] 

and  considering  the  relations  given  in  (13)  and  (16)  results  in 

(18)  PD(t+l)  =  *D(t+l,t)  [  Pj+  5Pj  ]  ®£(t+i,t)  +  r(t+ift)  Q(t) 

rT(t+i,t)  + 

+  E  |  6$>x(t/t-l)  x^(t/t-l)  6$  ] 

,  T  1 

+  E  |  6 A*  M'x( t / 1 -  1 )  xT(t/t-l)  M T6A'':  J  + 

[*  A  aT  T  1 

6A  M6x(t)  6x  (t)  M  6A  J 

where 

Pj  =  P(t)  -  P(t)  MT  [  MPMT-R(t)  j"1  MP(t) 

Round-off  errors  contribute  to  Pjj(t+1)  in  two  ways: 

1.  In  computing  the  estimator  x(t+l)  according  to  equation  (3) 
we  get  the  terms 

E  [  6<£>x  xT'6$T  ] 

E  |oA:'!  Mx  xT  MT£A',!  ] 

E  [dA*  M6x  6xT  MT6AV  ] 
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2. 


In  computing  the  covariance  matrix  P(t+1)  according  to 
the  Variance  Equation  (6)  we  get  the  terms 


6$  Pj6$  ^  and 


The  covariance  matrix  of  the  estimation  error  is  not  a  determinis¬ 
tic  function  of  the  time  anymore.  Its  elements  are  random  variables 
because  of  the  randomness  of  6$  and  6Pj.  Taking  the  expectation  (bias 
term),  we  obtain 


(19)  E  [  PD(t+l)  ]  =  P(t+1)  +  E  [  6$  Pj6$  T  ]  + 

E  [  6®x(t/t-l)  xT(t/t-l)  6$T  ]  + 

E  [  6A*  Mx(t)  xT(t)  MT6A*  ]  + 

T  1 

E  [  6A*  M6x(t,t-1)  6xT(t,t-l)  MT6A*  | 


where 


=  P(t)  -  P(t)  MJ 


MP(t)  M  +  R(t) 


-1 


MP(t) 


It  will  normally  be  sufficient  to  use  (20)  instead  of  (19): 

(20)  E  [  PD(t+l)  j  =  P(t+1)  +  e[s$  &(t/t-l)  xT(t/t-l)6$T  ] 


The  second  term  on  the  right  is  a  diagonal  matrix  with  the  elements: 

(21)  6P..(t+l)  =  £  e[»?.  ] 


6P„(t+  1)  =0  for  i  j  j 
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In  a  first  approximation  it  will  suffice  to  apply 


-K^]“  veK1- 

Thus , 

(22)  6P..(t+l)  ^  £x2eU.2  ] 

j  J  J 

One  Dimensional  Example 

The  system  is  determined  by 

x(t+l)  =  4  (t+ 1 ,  t)  x(t)  +  u(t) 

2(t)  =  x(t)  +  v(t) 

where 

2  2 

Eu  =  q  and  E  v  =  r 


The  expectation  of  the  variance  at  time  t+ 1  is  according  to  equation  (19) 
given  by 

E  [P(t+ 1)  =  E  |<J2(t+ 1 )  =  42(t+l,t)  + 

e[642  ]  +  q 

The  equilibrium  point  can  be  obtained  by  setting  p(t+l)  =  p(t)  =  p  : 

P  =  (iz  +  E642)  +  E[£2(t)  6;2  j 


4t)r 

p(t)+r  + 


'(t) 


or  (for  constant  t,): 

1 


S  ± 


f-  +  r  (q+  x2(t)  E6  42) 
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where 


S  =  r(;2  +  E6 1}  -  1)  +  q  +  x2(t)  E6£,2 

Existence  and  numerical  value  of  the  equilibrium  depends  obviously  on 
the  state  variable  x(t)  and,  therefore,  on  4- 

For  4  >  1  : 

lim  p  =  oo 
t  — *oo 

For  4=1  : 

lim  p  exists,  but  can  be  intolerably  large  if  the 
t  —*oo 

coordinate  system  cannot  be  chosen  properly. 

For  4  <  1  : 

lim  p  does  not  depend  on  the  state  variable  x(t). 
t  — *  oo 

If  4  is  a  function  of  time,  then  the  conditions  for  existence  of 
equilibrium  are  surprisingly  reduced:  4  can  be  larger  than  1,  if,  for 
instance,  4(1)  is  periodic  and  the  time  average  is  S  1,  But  it  is  not 
sufficient  to  require  only  that,  as  stated  in  reference  2 

0  <  6  ^  /4(t+l,t)  /  <  P,  <  oo 
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6.  MODIFIED  SYSTEM  MODEL 


Inspection  of  Equation  (3)  shows  that  propagation  of  injection  errors 
and  unperturbed  motion  of  the  vehicle  is  somewhat  intermixed.  The 
state  transition  matrix  determines  the  propagation  of  injection  errors 
as  well  as  the  unperturbed  motion  of  the  vehicle. 

If  we  define  two  transition  matrices,  one  for  the  perturbations  and 
one  for  the.  reference  trajectory,  then  the  estimation  problem  could  be 
separated  from  the  problem  of  integrating  the  equations  of  motion.  A 
system  model  is  shown  in  Figure  1: 

1.  £<J>(tj,t)  is  the  transition  matrix  for  the  perturbations  and 
determines  the  effect  of  injection  errors  at  an  arbitrary  time 
tj.  The  reference  trajectory  can  be  assumed  to  be  recompu¬ 
ted  every  time  our  knowledge  of  the  injection  conditions  im¬ 
proves. 

2,  M<J>(tj.t)  *s  t*ie  transit^on  matrix  for  the  state  vector  X  (t)  and 
determines  the  unperturbed  state  of  the  system  at  an  arbitrary 
time  tj.  The  elements  of  may  be  changed  every  time  our 
knowledge  of  the  actual  trajectory  improves  (the  perturbations 
tend  to  zero  if  the  estimation  procedure  converges). 

We  estimate  now  the  perturbations  x(t)  and  perform  the  computa¬ 
tion  of  the  state  vector  X  (t)  outside  of  the  estimation  loop.  We  assume 
that  knowledge  of  X  (t)  is  required  throughout  the  mission  to  decide  if 
corrective  maneuvers  should  be  made,  and  for  a  more  accurate  (not 
linearized)  representation  of  the  system  model. 

The  modified  error  equations  are  readily  found  with 

M<&D  = 

E^D  =  E^  +  ®  E^ 
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ESTIMATION  LOOP 


4 


4 


Figure  1.  A  trajectory  estimation  system  with  computation  noise; 
estimation  of  the  perturbation  vector  x(t.  ). 


[ 

[ 

[ 
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We  note  that  the  round-off  errors  which  are  generated  during 
numerical  integration  of  the  equations  of  motion  still  contribute  to  the 
estimation  error  xE(t): 

(23)  xD(t+  1)=  x(t  +  1/t)  +  8&(t+l) 
where 

8xi(t  +  1)  =  8E®ft(t/t-i)  +  8  A*  MxD(t/t- 1)  +  SM$x(t) 
or 

8$(t+l)  =  (SE$+  8M$)x(t)  +  SM$x(t)  +  8  A*  MxD(t/t  -  1) 


Repeated  application  of  the  recursion  formulas  above  leads  to: 

(24)  8$D(t)  =  ^  (  8E$  +  8m$  +  8A*M)I'1(8M$+8E$)x(t-i)  + 

i 

+  +  ®  A*  M)*  x  (t  -  i)  + 

i 

+  S  (8E$  +  SM$  +  8A*M)1‘1  SM$X(t-i) 

i 

and  shows  that  the  propagation  of  round-off  errors  from  previous  obser¬ 
vation  intervals  can  be  neglected,  if  8  ,  8E$  and  8A  are  sufficiently 

small.  In  a  first  approximation, 

(25)  8£D(t)  ~  (8m$+  Se$)  x(t-l)  +  (8m$+  8e$  +  8A*M)  x(t  - 1)  + 

+  8M*X(t-l) 

The  modified  covariance  matrix  is  given  by 
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E[PD(t+l)]=  E  ^E[xD(t+  1/t)  xD(t+  l/t)]J  = 

=  P(t+  i)  +  EfSg^Pjtt)  8e$t]  +  E[8M$x(t)xT(t)  8m#t  ]  + 

+  e[(8m$+  Se$)  £(t)  $T(t)  (8M$T  +  8e$T)  ]+  (26) 

+  e[8  A*M  xD(t)  xE(t)  MT  8a*T] 

Round-off  errors  in  the  filter  operation  A  *  (t)  can  safely  be  disregarded 
if  the  regular  estimation  error  x  (t)  converges  in  a  mean  square  sense, 
or ,  if 


lim  E  [7  x^J  = 
N— *oo 


We  write  therefore 

E[pD(t  +  1)]  =  P(t+  1)  +  E[8M$X(t)  xT  (t)  Sm$t  ]  + 

+  E  [(  8^  +  SE$)  x(t)  £T(t)  (8m$T  +  8  e$T)] 


(27) 


With  the  same  arguments  as  in  Section  5,  we  obtain  for  the  error 

in  element  P^  . 

D,  ii 


We  are  now  in  the  position  to  discuss  the  following  case:  Suppose 
one  would  like  to  reduce  the  effect  of  round-off  errors  by  estimating  the 
time  independent  injection  errors  x(tQ)  rather  than  the  current  pertur¬ 
bations  x(t).  The  state  transition  matrix  is  here 

E$(tj,  tQ)  =  1  and  8e$  =  0 
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But  in  computing  the  present  state  X(t)  each  time  a  new  measure¬ 
ment  is  made  or  each  time  a  decision  for  a  corrective  maneuver  is 
required,  the  equations  of  motion  have  to  be  integrated  from  t  to  t. 

The  variance  of  t  )  increases  significantly  because  of  the 

increased  integration  time  (Section  7)  and  may  finally  become  dominant. 
This  method  seems  favourable  only  if  ^$does  not  have  to  be  computed 
with  numerical  integration. 

If  we  are  forced  by  some  reason  to  determine  X(t^)  by  integrating 
the  equations  of  motion,  then  it  seems  advisable  to  modify  pur  estima¬ 
tion  procedure  so  that  (j^$  +  8^$)  becomes  part  of  the  estimation 
loop  (Figure  2). 

The  covariance  matrix  of  the  estimation  error  is  given  here  by 
the  diagonal  matrix: 

E[pD(t+  1)]  =  E  [8m$  X(t)  XT(t)  5  M®  ]  =  E  [fiXftj)  8XT(tl)] 

(29) 

of  0  0  0 

0  <7^  0  0 

0  0  of  0 

0  0  0  a  „2 

4 

The  of  are  equal  to  Var  fuj(t^)  j  and  given  in  Table  2. 
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igure  2.  A  trajectory  estimation  system  with  computation  noise; 
estimation  of  the  system  state  vector  X(t.  ). 


7.  EFFECT  OF  COMPUTATION  NOISE  ON  THE 
STATE  TRANSITION  MATRIX 


7.1  ERRONEOUS  PREDICTION  OF  THE  PERTURBATIONS 

A  method  for  computing  the  elements  of  the  state  transition 
matrix,  $  ,  has  been  described  in  (5).  The  elements  are  found 

jl* 

there  with  six-fold  numerical  integration  of  the  perturbation  equa¬ 
tions 


x  (t)  =  F  (t)  x  (t) 

under  the  appropriate  initial  conditions.  We  choose  this  method  as  an 
example  for  studying  the  round-off  error  propagation,  because  it  is 
quite  generally  applicable  and  requires  especially  high  computation 
accuracy  (worst  case  example). 

We  are  not  interested  here  in  a  complete  round-off  error  analysis 

where 

the  method  of  integration  used, 

the  form  of  the  equation  to  be  integrated,  and 

the  program  organization 

would  have  to  be  considered  in  detail.  We  restrict  our  study  to  a  sym¬ 
metric  round-off,  Heun's  method  of  integration,  and  assume  optimal 
scaling.  We  also  do  not  analyze  the  effect  of  the  equation  form  on  error 
propagation.  However,  we  keep  in  mind  that  canonical  transformations 
(separation  of  the  variables)  could  probably  always  be  used  to  reduce 
the  accumulated  round-off  error.  Each  equation  of  the  system  could  be 

integrated  separately  by  means  of  the  trapezoidal  rule  and  the  accumu- 

1  ^M 

lated  round-off  error  could  be  made  proportional  to  or  — for 

each  variable. 

We  intend  to  derive  relations  between  the  statistical  moments  of 
the  round-off  errors  in  the  state  transition  matrix  and  computation  par¬ 
ameters  like  wordlength  and  number  of  integration  steps. 

It  seems  desirable  for  our  purpose  to  obtain  mean  and  variance 
in  closed  form  rather  than  by  simulation  methods,  A  circular  orbit 
with  the  state  variables  r,  0  ,  r,  9  was  chosen,  therefore,  as  an  example 

The  error  Equations  (33)  and  (34)  can  be  solved  in  closed  form. 
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Generalization  to  non  circular  orbits  and  6  state  variables  can  readily 
be  accomplished  (Appendix  A)  but  not  in  closed  form.  We  would  have  to 
find  polynomial  approximations  for  the  solution  of  the  adjoint  system 
(34).  However,  the  results  in  this  section  can  also  be  used  to  approxi¬ 
mate  the  round-off  error  in  the  general  case. 

The  error  variables  are  defined  as  the  deviations  of  the  actual 
from  the  reference  trajectory: 


=  r  -  r 


R 


=  6-0 


R 


where  r  and  0  are  the  polar  coordinates  of  the  vehicle  with  the  origin  in 
the  center  of  attraction. 

From  the  equations  of  motion 


r 


K 

'  “T 


+  r  0 


2 


0  = 


(30) 


we  find  readily  the  perturbation  equations  (perturbations  are  assumed 
to  be  caused  by  injection  errors  only)  with  Taylor  expansion  of  the  right 
side  of  (30).  The  perturbation  equations  are 

xl  (t)  =  x3  (t) 

x2  (t)  =  x4  (t) 

x^  (t)  =  3co  ^  Xj  (t)  +  2  r  w  (t)  ;  w  =  9  =  constant  (31) 

x4  (t)  =  -  -  x3  (t) 

and  could  be  integrated  in  closed  form.  But  we  integrate  numerically 
with  Heun's  method*  [4]  to  study  the  propagation  of  round-off  errors. 
We  compute  first  the  auxiliary  values 

*Heun's  method  is  less  accurate  than,  for  instance,  the  Runge 
Kutta  method,  but  simpler  to  analyze.  Its  accuracy  seems  sufficient 
for  our  objective  "to  study  the  effect  of  wordlength  on  the  accuracy  of 
the  state  transition  matrix." 
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xiS<‘j)  =  Xi<tj-l)  +  At'i  .  Xn,*j-1)’ 


and  then 


At 


xi(tj)=xi(tj_l)+“2-fi 


i  x  (t.  x  (t.  ) i  +  f.  r  x  "'(t.), ....  x 

L  1  J-l  n  j- 1  J  i  l  1  '  J  n v  j' 


i  =  1,  .  .  .  ,  n  and  j-l,  ,  N 


The  determination  of  round-off  errors  can  be  considerably  simplified, 
if  we  assume  that  At  is  a  small  quantity  so  that  multiplication  by  At 
shifts  the  doubtful  figures  in  A  (x^  .  .  .  x^)  into  those  which  are  dropped 
in  the  last  round-off  step.  We  have  the  identity 


t  • 


x 

n 


where  the  bars  indicate  round-off  values. 

If  we  are  not  interested  in  higher  powers  of  At,  then  we  may  write 
for  the  jth  step: 


Llj 


j-l 


it  x 


3  j-l 


+  2  e. 


x_  .  =  x_ .  +  At  x .  .  +  2  e 

2  j  2  j - 1  4  j-l  2 


x,  .  -  x,  .  .  +  At  [  3w  2  x  +  2  rw  x  ]  +  2  e 
3  J  3 j-l  \  1 j-l  4 j-l  I  3 


(32) 


x.  .  =  x A  .  ,  -  —  x0  .  +  Z  eA 

4j  4  j- 1  r  3  j- 1  4 


which  we  have  to  compare  with 


x  .  =  x,  .  ,  +  At  x-  . 
lj  lj-1  3  j-l 


4 j  "  ^4  j-l 


3  j-l 
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Introducing  the  error  variables 


u.  .  =  x.  .  -  x.  . 
iJ  iJ 


we  can  write  the  error  equation  for  step  j 

u,  .  =  u. .  +  At  [  x.  .  -  x  .  1  -  2  e 

lj  lj-1  L  3  j- 1  3  j- 1  J  1 

u .  ■  -  u .  .  +  At  f  -  (— )  x-  .  .  +  (— )  x  1-2 

4j  4 j-1  L  r  3 j-1  r  3 j-1  J 


or 


u.  .  =  u.  .  +  At  u  .  -  2  e 

lj  lj-1  3  j- 1  1 


u_.  .  =  u_  .  +  At  u  .  .  -  2  e, 

2  j  2  j - 1  4  j-1  2 


u,  ..  =  u,  .  +  At  3  w  u.  .  +  2  r  oj  u  -  2  e 

3 J  3 j-1  l  1 j-1  4 j-1  J  3 


u.  .  =  u.  .  -  At  —  u  j  -  2  e 

4  J  4  j- 1  l  r  3  j - 1  /  4 


To  solve  Equation  System  (33)  we  introduce  the  adjoint  system 

X,.  =  X,  .  -  At  •  3w  2  X 

lj  lj-1  3  j 


(33) 


2j  = 

X_  . 

2  J-1 

X.  . 

-  At 

f  xi  • 

3  j  = 

3  j-1 

l  ij 

X  .  . 

-  At 

(  X 

4  j  = 

4  j-1 

(34) 


3  j 
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•V 

f 


and  infer 


(X  1  j  U1  j  -  X  1  j-l  ulj-l)  +  (x2ju2j  '  x2j-lu2j-l)  +  (x3ju3j"  X  3  j  —  1  ^j-l*  + 


+  (x4ju4j'  x4j-lu4j-l)  "  ■  2(Xljel+  X2je2  +  X3je3  +  X4j*4* 


Summing  over  j  from  1  to  N  we  obtain 


X1  U1  +  X  2  U2^W  +  X3  U3  +  X  4  U4 


2Z<X1J*1  +  X2j,2+X3j*3+X4j%)  <35> 

j  =  l 


where 


t,  .  =  NAt  =  total  time  of  integration 


At  =  integration  step 

Treating  the  e's  as  independently  and  uniformly  distributed,  we 
can  write  down  the  variance  of  (35) 

4  1  N 

V“  2  Xi(,M>  =42(X  lj  +  X2j+X3j  +  X4j)t’2  <36> 

L i=l  J  j=l 

where  <r  ^  designates  the  round-off  error  variance  in  one  integration  step. 
The  sums  on  the  right  can  be  looked  upon  as  Rieman  sums  and  can  be 
replaced  by  the  integrals 

rJL  i  ?  rtM  ,  _ 


->  i  2  r 

Var  Z  xi(tM>  ui(tM>  =¥rl 

-  i=l 


X^t)  +  V  22(t)  +  X  32(t)  +  X  42(t)J  dt 


order  to  get  the  variances  in  u^  separately,  one  has  to  impose  on 


X  .  (t)  the  terminal  conditions: 
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1. 


w  =  ‘  •  kw  =  w  =  vy=° 

2-  X2^tM)  =  1  ’  =  X3^M^  =  X4^tM)  =  0 


3'  X3^tM)  1  ’  =  X "  VW  “  0 


(38) 


4-  ~  1  ’  S^M*  ”  X2^M^  ~  X3^M^  "  0 

SOLUTION  OF  THE  ADJOINT  SYSTEM 

We  can  write  the  adjoint  system  in  differential  form: 


Xj  =-3  «‘X3 


X2  =  ° 


x3  =  -  Xi  \4 


(39) 


X  =  -  X  -  2  r  «  X 
4  2  3 


and  obtain 


X  =-C  -  —  C.  t  -  3wC  e  wt  +  3WC,  e_Wt 

1  r  4  r  1  2  3 


X2  C1 


1  lAi  t  -  {J  t 

X3=7 i-*C2e  +  C3 


(40) 


X  =  C  .  -3C,t-2rC_ewt  +  2rC,  e‘Wt 
4  4  1  2  3 


For  the  four  sets  of  terminal  conditions  we  get  four  sets  of  solutions 
(Table  1). 


25 


Table  1.  Solutions  of  the  adjoint  system  for 
various  terminal  conditions. 
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The  resulting  variances  of  the  round-off  errors  at  an  arbitrary 
time  tj^j  are  tabulated  in  Table  2.  In  a  first  approximation,  the  vari¬ 
ances  in  the  coordinates  x^ ,  X£  and  x^  increase  proportional  to  the 

number  of  integration  steps  N,  but  the  variance  in  x0  increases  pro- 

3  i 

portional  to  N  . 


The  error  in  the  state  transition  matrix  $  (t  ,  t)  now  can  be 

£  1 

characterized  by 


[*£«]- 


and 


(41) 


where 


equal  to 


Var[ui(tM>] 


and  given  in  Table  2. 
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7 .2  ERRONEOUS  PREDICTION  OF  THE  STATE  VECTOR 

Restricting  ourselves  again  to  an  in-plane  trajectory  of  a  satellite, 
w.e  have  to  integrate  two  second  order  differential  equations  of  the  form 
(30). 

Setting 

r  =  Xj  ,  0  =  X2  ,  r  =  X3  ,  0  =  X4  , 


we  can  write  the  equations  of  motion  as  a  system  of  first  order  equations 


X,  =  X, 


X2  =  X4 


*3  *  -  +  X1  x4  -  fl  <xr  x2,  X3,  X4) 

X1 


X3X4 


=  f2  (Xj ,  X2  ,  x3  ,  x4) 


(42) 


and  obtain  the  error  equations  similar  to  those  in  Section  7.  1 


u  .  .  =  u.,  .  ,  +  At  u  ,  .  ,  2  e 


lj  “  lj  -  1 


3j  -  1 


u  =  u 
2j  2j  -  1 


+  u  4j  -  1  '  2  e2 


U  3J  =  u 3j  -  1  +  At  [(^T  +  X4  )  “  1  +  '2X1  X4>  "4_ 


2  e. 


“4j  =  u4j  -  1  +At[( 


hl±ri 

xf  1  1 


2  €  . 


(43) 
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The  only  difference  is  that  the  coefficients  are  functions  of  time 
here  and  not  constants.  This  does  not  complicate  matters  too  much 
because  coefficient  accuracy  is  not  critical.  Nominal  values  of  the  tra¬ 
jectory  suffice.  We  may  go  even  one  step  farther  and  restrict  ourselves 
to  an  approximately  circular  orbit  and  determine  the  coefficients  by 
setting 


=  r  -S  constant 

X3  »  0 

2 

X.  =  ui  S  constant 
4 


The  equation  system  (43)  reduces  to  (33)  and  can  be  solved  with  the 
elementary  methods  described  in  Section  7.1. 

The  error  vector 


8  X  (tj)  =  8M$8x(t) 


now  can  be  characterized  by 

E  [8X(t1)]  =  0 


and 


E  [sx(t1)  8  XT  (t  j  )] 


0 

0 


0 

0 


(44) 


where  the  a are  equal  to  Var  £u  ^  (t^)  J  an<^  given  in  Table  2. 
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8.  EXAMPLE:  SATELLITE  TRAJECTORY  ESTIMATION 
BY  THE  ONBOARD  DIGITAL  COMPUTER 


We  assume  an  approximately  circular  orbit  in  the  terrestrial  field. 

The  actual  in-plane  trajectory  is  perturbed  around  a  precomputed  ref- 

7 

erence  trajectory  (r  =  2.2  •  10  feet).  Observations  are  made  period¬ 

ically  and  used  in  the  spaceborne  digital  computer  to  estimate  the  current 
position  and  velocity  in  the  trajectory  plane  according  to  Figure  2. 

We  are  interested  in  finding  the  estimation  error  caused  by  the 
finite  wordlength  in  the  digital  computer,  and  to  separate  the  effect  of 
wordlength  K , 

observation  interval  tw,  and 

number  of  integration  steps,  N,  in  one  observation  interval 
on  the  error  distribution.  We  know,  that  the  error  is  mainly  determined 
by  the  round-off  errors  in  the  numerical  integration  of  the  equations  of 
motion.  The  variances  of  the  accumulated  error  are  given  in  general 
form  in  Table  2.  Figure  3  shows  the  error  propagation  in  the  4  co¬ 
ordinates  r,  9,  r,  0  for  our  example  as  a  function  of  the  integration 
time  tj^  and  for  a  0.  5  seconds  integration  interval. 

The  error  grows  quite  rapidly  with  N(or  t^  =  NAt).  The  result  is 
not  in  agreement  with  the  commonly  ^  assumed  error  growth  propor¬ 
tional  to  N.  It  does  not  seem  realistic  to  assume  that  the  round-off 
errors  grow  proportional  to  N  independently  of  equation  structure,  inte¬ 
gration  method,  coordinate  system,  etc.  It  is  noted  in  Table  2  that  the 
error  grows  proportional  to  N  (in  agreement  with  ref.  6),  but  only  for 
three  of  the  four  state  variables  and  for  a  very  short  integration 
time. 

The  curves  in  Figure  3  are  independent  of  the  computation  accu- 
2  2, 

racy.  If  we  multiply  o,  with  the  variance,  a  ,  of  the  round-off  error 

1  2 
in  one  integration  step,  we  find  the  actual  error  variance  o.  -  Table  3. 

We  look  now  for  a  geometrical  representation  of  our  error  co- 
variance  matrix,  and  introduce  the  ellipse  of  concentration.  The  ellipse 
(Figure  4)  is  defined  as  a  curve  enclosing  all  points  which  deviate  from 
the  reference  point  with  a  probability  at  most  equal  to  P. 
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Figure  3.  Accumulation  of  round-off  errors  during  the  numerical 
integration  of  a  system  of  first  order  differential 
equations.  Integration  step  0.  5s. 
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Table  3.  Variance  of  the  accumulated  round-off  errors  as  a  function  of  wordlength  and  integration  time 


Figure  4.  Ellipse  of  concentration  as  a  performance  measure  for 
the  digital  computer  in  midcourse  navigation. 


Let  x  and  y  be  cartesian  coordinates  with  respect  to  the  centre 
of  the  ellipse  (on  the  reference  trajectory) 


x  =  Ar 
y  =  rA  6 

where  Ar  and  A0  are  small, 
then  the  equation  of  the  homothetic  ellipses  is 


2 

x 


2  2 
r  °  2 


2c 


(45) 


2  2  2 
where  a,  and  <7,  are  given  in  Table  2  and  Table  3.  The  constant  c  is 

2  *  ^ 

X  distributed  with  2  degrees  of  freedom.  The  probability  of  falling  out¬ 
side  of  the  ellipse  is  obviously  given  by 


P  =  P(X2  >  2c2). 


(46) 


It  is  noted  that  in  our  example  (Figure  4)  the  major  axes  of  the  ellipse 
always  lie  in  the  direction  of  r  and  perpendicular  to  r  (downrange).  But 
a  coordinate  transformation  in  equation  system  (42)  would  change  length 
and  direction  of  the  axes. 

The  effect  of  wordlength  on  the  error  ellipse  is  shown  in  Table  4. 

The  area  covered  by  the  ellipse  grows  in  our  particular  example  from 

2  2  ..  io  2 

3.  3  •  10  feet  to  the  discouraging  magnitude  of  1 .  7  •  10  feet  if 

we  require  a  probability  of  99-9  percent  of  not  falling  outside  the 

ellipse.  Obviously,  we  may  generate  any  size  ellipse  and  require  any 

wordlength  for  the  digital  computer,  if  we  can  persuade  ourselves  to 

more  or  less  optimistic  probability  figures  (Figure  5). 
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r 

I. 


r 


K 

a  [feet] 

b  [feet] 

F  [  feet  ] 

24 

1.8  x  104 

3.05  x  105 

1.7  x  1010 

28 

1.12  x  103 

1.91  x  104 

2.1  x  107 

30 

0.282  x  103 

0.47  x  104 

1.35  x  106 

32 

0.7  x  102 

1.19  x  103 

8.4  x  104 

34 

1.76  x  101 

2.98  x  102 

5.3  x  103 

36 

4.4 

7.45  x  101 

3.3  x  102 

Table  4.  Half  axis  and  area  of  the  99.9  percent  error 
ellipse  K  =  word  length,  observation  inter¬ 
val  10  minutes,  approximately  circular  orbit 
with  2.2  x  10?  feet  radius. 


36 


ir  ellipses  for  varies  error  p 
.gits,  observation  interval  10 
liar  orbit  with  2.  2  x  10^  feet 


9.  SUMMARY 


The  error  ellipse  (ellipse  of  concentration)  has  been  introduced 
as  a  performance  measure,  and  its  axes  and  area  have  been  expressed 
in  terms  of  the  computer  wordlength,  the  time  interval  between  obser¬ 
vations,  and  the  number  of  integration  steps  for  a  particular  class  of 
estimation  systems.  Characteristic  of  the  system  is  that  the  equations 
of  motion  are  solved  by  numerical  integration  and  not  analytically.  The 
conclusions  about  the  growth  of  the  estimation  error  are,  therefore, 
rather  pessimistic.  A  computer  designed  to  operate  satisfactorily  in 
our  example  certainly  would  also  be  appropriate  for  more  elegant  navi¬ 
gation  systems . 

As  a  sideline  to  our  main  problem,  it  was  necessary  to  obtain 
some  information  about  the  round-off  error  propagation  in  numerical 
integration  of  first  order  differential  equations  (second  order  equations 
can  readily  be  reduced  to  first  order  equations).  The  error  moments 
have  been  found  by  means  of  a  first  order  perturbation  method  and  have 
been  expressed  in  terms  of  the  adjoint  solutions  of  the  error  equations. 

In  general,  the  error  variances  grow  with  some  power  of  N.  The  power 
depends  on  the  integration  time,  number  of  variables  in  the  equation 
system,  and  on  the  form  of  the  equations  (coordinate  system).  The  co¬ 
efficients  of  the  differential  equations  have  been  made  approximately 
constant.  Otherwise,  the  adjoint  equations  cannot  be  integrated  in  closed 
form. 

If  we  are  satisfied  with  first  approximations  of  the  estimation  error 
and  if  we  may  assume  convergence  of  the  regular  estimation  error  to  a 
moderately  small  equilibrium  (steady  state),  then  the  problem  of  finding 
a  performance  measure  for  the  digital  computer  during  midcourse  navi¬ 
gation  is  reduced  to  the  task  of  determining  the  round-off  error  propa¬ 
gation  during  numerical  integration  of  the  equations  of  motion  and 
perturbation.  Table  2  gives  closed  form  expressions  for  the  accumulated 
round-off  error  variances  assuming  approximately  circular  orbit  and 
in-plane  two  body  motion  (Section  7.  1  and  7.2). 
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The  effect  of  limited  computer  accuracy  on  the  estimation  error 
can  be  determined  with  the  Formulas  (28)  and  (29).  Both  formulas  are 
approximations.  Neglected  are: 

a.  Propagation  of  round-off  errors  from  previous  observation 

intervals  against  the  error  contributions  of  the  last  interval. 

Round-off  errors  from  previous  observation  intervals  would 

contribute  only  in  pathological  cases  where  X(t.)  decreases  so 

2  1 

rapidly  with  T  that  26<j>  X(t-2)  is  not  neglegible  against 
6  4»X(  t  - 1) . 

b.  The  effect  of  the  erroneous  filtering  operation  A  ‘  (t)  +  SA  *. 

It  can  safely  be  done  only,  if  the  regular  estimation  error 
"x(t)  converges  with  the  increasing  number  of  measurements: 
lim  E  lx"  St"1 1=0 

N— *oo 

The  errors  in  the  state  transition  matrices  „4or  are  associ- 

M  E 

ated  in  (28)  and  (29)  with  the  perturbation  vector  x(t)  or  with  the  state 
vector  X(t).  Even  very  small  6<j>may  contribute  significantly  to  the 
estimation  error ,  if  x(t)  or  X(t)  are  sufficiently  large.  A  typical  ex- 

12 

ample  is  the  transition  matrix  of  injection  errors  for  a  circular  orbit  . 
The  matrix  contains  elements  which  grow  with  t.  Unbounded  perturba¬ 
tion  components  x^(t)  and  unbounded  solutions  of  the  variance  equation 
must  be  expected. 

We  notice,  that  computation  noise  may  increase  the  uncertainty 
in  our  knowledge  of  the  system  state  ,  even  if  a  large  number  of  observa¬ 
tions  have  been  made.  We  found  that  the  covariance  matrix  of  the  com¬ 
putation  noise  appears  to  be  additive  to  the  covariance  matrix  of  the 
random  force  u(t)  in  (18).  We  may  define,  therefore,  a  generalized  co- 
variance  matrix 

Q '(  t)  =  Q(t)  +  [pD(t)-P(t)j 

and  recognize  that  Q  (t)  may  be  unbound.  We  conclude  that  theorem  4 
in  reference  2  does  not  apply  for  dynamic  systems  which  are  corrupted 
by  computation  noise. 
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APPENDIX  A 


GENERAL  ADJOINT  DIFFERENTIAL  EQUATIONS 
FOR  THE  ROUND-OFF  ERROR 


The  equations  of  motion  are  generally  given  by  a  nonlinear  system 
of  second  order  differential  equations.  We  write  it  as  a  system  of  first 
order  differential  equations: 


V  x4 

x2=  x5 

V  X6 


4  = 

fl 

<X1  , 

x2  . 

....  x6> 

5  = 

f2 

(Xi  ■ 

x2  - 

x6) 

6  = 

f3 

<X!  ■ 

■  X2  ’ 

....  x6) 

(47) 


and  linearize  by  Taylor  series  expansion  around  the  reference  trajectory. 
Introducing  the  perturbation  coordinates 

Xi=  Xi  -  Xi,  R  •  1  =  1 . 6 

we  obtain  for  the  perturbation  equations: 

*1=  *4 

x2=  x5 

*3  =  x6 

6 

Vi 

i=  1 


(48) 


dx{ 


li 


(x.) 


i  =  1 
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•  -  v  fli 
"5  =  A  3Xi 


i = i 
6 


*e-Z 


i  =  1 


x. 


df3 

dXi  xi 


X  F2i(Xi' 


i  =  1 
6 


=  £  F3i(x.) 


i  =  1 


Numerical  integration  on  a  digital  computer  using  Heun's  method  leads 
to  round-off  errors  which  can  be  characterized  in  step  j  by  the  error 
equations: 

ulj  =  ulj-l  +  At  u4j-l  '  2el 
U2j  =  u2j  -  1  +  AtU5j-l  -  2e2 


U ,  •  =  u  ,  •  ,  +  At  u  /  .  ,  -  2e 


3j  3  j  -  1 


u4j  =  u4j  -  1  +  At 


%=  "5j-l  +  At 


6j  -  1 

6 


I 

i  =  1 

6 

X 

i  =  1 


.3 


dFH 

a  xi  U  ij  -  1  2  e4 

*F2i  u 

dxi  ij'1  65 


(49) 


U6j  =  %  -  1  +  At 


V5  ^F3i 

>  - —  u.  -2  e, 

2ij  a  xi  1  6 


i  =  1 


where  the  e.  designate  the  error  in  one  integration  step  and  u. .  = 
Xi(t.)  -  x.(t.). 
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We  solve  the  Equation  system  (49)  by  means  of  the  adjoint 
equations 


xu  =  N-i  -At  E 


aF 


K1 


K  =  1 


axi  (K  +  3),  j 


X2j  '  V. 


•J 

1 
K  =  1 


aF 


K2 


dx2  (K  +  3),  j 


X3j  -  1  - 


X 

K  =  1 


aF 


K3 


dx3  (K  +  3),  j 


(50) 


X4j  =  X4j-1  -  At<  X-  + 


X3J  = 


X5j  -  1  "  At  '  + 


V  aFK4 

Z.  ax4 


K  =  1 
3 


4  (K  +  3)  j 


V  K5  \ 

Zj  ax=  (k+3)j 

K  =  1  3 


^  ~  Xv  ■  ,  -  A  t  \  X.,  + 


6j  “  *6j  -  1 


3j 


y  _?Ik6  x  ) 

A  ax6  (K  +  3)  j ^ 


K  =  1 


and  obtain  the  variance  of  the  round-off  error  in  the  coordinates 
5V<W  at  time  readily  with  the  argumentations  in  Section  7.1: 


r-  6 


Var 


2  XiV>M(tM> 


i  =  1 


4  or 


2  6 


At 


/  X  Xf(t)dt-  <51> 


o  i  =  1 
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